Interplay of inter-chain interactions and exchange anisotropy: 
Stability of multipolar states in quasi-lD quantum helimagnets 



S. Nishimoto, 1 S.-L. Drechsler. l! [j R. Kuzian, 1,2 J. Richter, 3 and Jeroen van den Brink 1 

1 IFW-Dresden, P.O. Box 270116, D-01171 Dresden, Germany 
'^Institute for Problems of Materials Science NASU, Krzhizhanovskogo 3, 03180 Kiev, Ukraine 
3 Universitat Magdeburg, Institut fur Theoretische Physik, Germany 
(Dated: March 11, 2013) 

We quantify the instability towards the formation of multipolar states in coupled spin- 1/2 chain 
systems with a frustrating J1-J2 exchange, in parameter regimes that are of directly relevance 
to edge-shared cuprate spin-chain compounds. Three representative types of inter-chain coupling 
and the presence of uniaxial exchange anisotropy are considered. The magnetic phase diagrams are 
determined by Density Matrix Renormalization Group calculations and completed by exact analytic 
results for the nematic and dipolar phases. We establish that the residual couplings strongly affect 
the pitch of spiral states and their instability to multipolar phases. Our theoretical results bring to 
the fore novel candidate materials close to quantum nematic/triatic ordering. 



In a system with frustrated magnetic interactions en- 
tirely new ground states (GS) can emerge from the en- 
suing competition. The geometric frustration of classical 
Ising spins on a pyrochlore lattice, for instance, results 
in the famous spin-ice state, the excitations of which are 
magnetic monopoles [1 . In frustrated quantum magnets 
equally exotic states such as spin liquids, valence-bond 
crystals or nematic phases, can occur [2J. In quantum 
spin chain systems, in particular, the competition be- 
tween short and longer-range magnetic couplings is a 
common source of frustration, a canonical example of 
which is the Ji-</ 2 spin-1/2 chain [2J. Having antiferro- 
magnetic (AFM) next nearest-neighbor (NNN) interac- 
tions ( J 2 > 0), it is frustrated for any sign of the nearest- 
neighbor (NN) coupling (Ji). In the classical J\-J% spin 
chain the competing interactions generate a helimagnetic 
state but in a single chain quantum fluctuations destroy 
the long-range helical order for any value of J\. For suf- 
ficiently high magnetic field, for any value of J\, the FM 
state takes over and the system's magnons, its propa- 
gating spin-flips, become its exact single-particle excita- 
tions. The exchange parameters Ji and J 2 determine 
the magnon dispersion and, in particular, the interaction 
between them. An AFM interaction leads to a repulsion 
between magnons, whereas a FM interaction results in an 
attraction, which favors the formation of magnon bound 
states. For a frustration ratio a=— J 2 /Ji > 0.367 an in- 
teresting and intensely studied nematic state can occur, 
which may be thought of as a condensate of 2-magnon 
bound states [3"HT5] characterized by a quadrupole spin 
order with a non-zero anomalous average (S^Sj~). For 
1/4 < a < 0.367 also 3-, 4- and even higher magnon 
bound states can condense, resulting in a rich phase dia- 
gram with quite a number of exotic magnetic multipolar 
phases (MPPs). 

These theoretical developments have stimulated an 
experimental quest to find multipolar condensates in 
quasi one-dimensional (ID) magnetic materials, in par- 
ticular in spin s = 1/2 systems consisting of edge-sharing 



copper-oxide chains, such as LiVCu04 (in cuprate nota- 
tion =LiCuV04 in traditional chemical notation)) |10l 
EME Li 2 ZrCu0 4 HZ1 US], Ca 2 Y 2 Cu 5 Oio [H ED], 
PbCuS0 4 (OH) 2 (HI H2], Rb 2 Cu 2 Mo 3 Oi2 [S3] and 
Li 2 Cu0 2 |24] [25] . In these systems J\ is intrinsically 
FM and J 2 can be of comparable strength, but AFM. In 
real 3D materials, however, a magnetic inter-chain (IC) 
interaction is unavoidably present. Due to the fragility of 
purely ID bound-states, AFM IC interactions can pose 
a very relevant perturbation to a multipolar state, even 
when the coupling strength is (very) small [55]. To es- 
tablish the consequences of this key ingredient for the 
stability of MPPs we consider here the three most com- 
mon types of IC couplings J IC that are encountered in 
the quasi-lD edge-shared cuprates mentioned above (one 
perpendicular IC coupling and two different types of skew 
ones, see Fig. [lj and determine the boundaries of the 
magnetic phase diagram numerically by Density Matrix 
Renormalization Group (DMRG) calculations and ana- 
lytically by hard-core boson (HCB) |25ff2l)] ) and spin- 
wave (SW) [20] approaches. On top of this we consider 
also the presence of a uniaxial exchange anisotropy A — 1 
for the NN coupling along the chains, which is the leading 
anisotropy term in edge-shared chain cuprates |30If3"3"] . 
We show that the stability of MPPs is strongly affected 
by the strength of the AFM IC couplings and depends 
on the precise type (geometry) of this coupling, which 
may also largely affect the pitch of the spiral state. A 
small easy-axis exchange anisotropy, however, enhances 
the stability of MPPs dramatically, also in the presence 
of IC coupling, since it enhances the attraction between 
magnons. From the material's viewpoint, our theoret- 
ical results bring to the fore linarite, PbCuS04(OH) 2 , 
as a promising candidate compound with a triatic MPP, 
which can be stabilized by its sizable exchange anisotropy 
and confirm the closeness of LiVCu04 to quantum ne- 
maticity. 

The relevant Hamiltonian H = Hid + Hie encom- 
passes the frustrating magnetic interactions along the ID 
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Figure 1: (a) Competing NN and NNN exchange Ji and J2, 
respectively, along a chain. Coupling between different chains: 
(b) perpendicular coupling Jq C (e.g., LiVCuCU), (c) skew (di- 
agonal) coupling j( c (e.g., PbCuS04(OH)2) and (d) skew 
NNN coupling between shifted chains Jl c (e.g., Li2Cu02). 
The effect of J IC is considered in both 2D and 3D. 



chain in the presence of an external magnetic field h and 
a small uniaxial exchange anisotropy A — 1 

Hid = J]] [— Sn,i • S„^ +i + aS n> i ■ S ni j+a] 

n,i 

- [( A - + hS* ti ] , (1) 



where n labels the chain and i the position of the spins 
along the chain. Neighboring chains n and m interact via 
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where r = corresponds to a perpendicular IC coupling 
and r = 1, 2 refer to skew IC couplings, see Fig. [T] We 
use I Ji I as the energy unit of all coupling constants in H. 

To determine the nature of the magnetic GS and its 
dependence on the frustration a, the different types of 
IC exchange J IC and the exchange anisotropy A — 1, we 
employed the DMRG method pH] with periodic bound- 
ary conditions (PBC) for all directions. This method 
is not restricted to purely ID and can also be used for 
2D [55J [35] and 3D [551 [59] systems, although the sys- 
tem size is limited, e.g., up to about y/~N x y/N x L = 
y/lQ x y/lO x 50 for spin Hamiltonians. We kept p w 
800 — 5000 density-matrix eigenstates in the renormal- 
ization procedure. About 100 — 300 sweeps are nec- 
essary to obtain the GS energy within a convergence 
of 10~ 7 Ji for each p value. All calculated quantities 
were extrapolated to p — > 00 and the maximum error 
in the GS energy is estimated as AE/Ji ~ 10~ 4 , while 
the discarded weight is less than 1 x 10~ 6 . Under the 
PBC, a uniform distribution of (Sf) may give an indica- 
tion to examine the accuracy of DMRG calculations for 



states [5 t z ot > (NL - 10) /2)\ the GS energy can be ob- 
tained with an accuracy of AE/Ji < 10~ 12 by carrying 
out several thousands sweeps even with p sa 100 — 800. 

We considered systems with different lengths: L = 
16 - 64 (24 - 96) for 3D (2D) and adopted power laws to 
perform a finite-size-scaling analysis. From this we ob- 
tained the saturation field h s in the thermodynamic limit 
L — > 00. As a result, we obtain h s with high accuracy. 
In addition to DMRG we have also applied an analytic 
HCB-approach and the linear SW approach [261 |2"T] to 
provide exact results for the nematic and dipolar phases. 
In addition, some of the calculated magnetization curves 
have been cross-checked by exact diagonalization. 

The simplest case, relevant for, e.g., L1VCUO4 and 
Li(Na)Cu202, is the situation of unshifted neighboring 
chains and a perpendicular inter-chain exchange Jg C , see 
Fig.[T^. In this case spirals on NN chains are only weakly 
affected by an AFM IC coupling [37 - on a classical level 
the pitch of the incommensurate (INC) spiral state is not 
affected by J^ . This is in stark contrast to the effect 
of skew AFM J\ G and j\ c , which can strongly reduce 
the pitch. A typical magnetization curve for a=0.5 and 
A = 1, for a nematic phase, is shown in Fig.[2l The height 
of the magnetization steps AS* Z =2 when J^ c /a=0.1, is 
the direct signature for 2-magnon bound states. A larger 
value of the IC coupling suppressed these bound states, 
as is clear from the magnetization curve for Jq /cv=0.2 
where the steps correspond to AS' Z =1. So in the isotropic 
case, where A = 1, a rather weak critical IC of a few 
percent destroys the nematic phase in favor of the usual 
conical ordering. The critical value for Jq C /(a = 0.5) 
amounts 0.188/0.088 in 2D/3D, respectively. The full 
phase diagram [3H] is shown in Fig. [3] where the phase 
boundaries are extracted from the kinks in the calculated 
saturation field h s as a function of Jq C , as shown in Fig. 
|3ja-c). Clearly, the 3- , 4- , and higher multimagnon 
MPPs are even stronger affected by the IC interaction. 

Allowing for a finite uniaxial exchange anisotropy 
A — 1, the leading-order anisotropy that is of immedi- 
ate relevance to quasi- ID cuprates |30[ I3T] affects the 




spin systems. Typically, (S z ) 
1 x 10~ 3 in our calculations. 



- S% ot /(NL) is less than 
Note that for high-spin 
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Figure 2: Magnetization vs. magnetic field for a 2D arrange- 
ment of four chains with iV = 24 sites each, with a perpen- 
dicular IC coupling Jl c (cf. Fig. []}>), a = 1/2 and A = 1. 
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stability of the MPP substantially. Fig. [3] shows that for 
a = 1/2 an anisotropy A— 1 of just 0.1 increases the criti- 
cal IC coupling by a factor of ~1.6, and thus significantly 
enhances their stability region. 

Our analytical approach to calculate the phase bound- 
ary between the 1- and 2-magnon instabilities relies on 
first deriving the saturation fields of these two instabil- 
ities: h Sl i and h S 2 respectively. Requiring them to be 
equal then renders the equation for the critical IC cou- 
pling as a function of anisotropy and frustration param- 
eters. The saturation field h Sj i of the INC phase on the 
1-magnon side is exact already within SW theory: 



(4a - l) 2 
8a 



IC , 7 IC 



G/o° + l</o°l)-(A-l), (3) 




where N\c denotes the number of IC neighbors (i.e. for 
Jg C in 3D and 2D, Nic = 4 and 2, respectively). In 
the Supplementary Material this expression has been fur- 
ther generalized to include next NN and IC exchange 
anisotropics [57]. Therein we have shown also that the 
critical value of Jq C of perpendicular IC depends only 
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Figure 3: (a-c) Saturation field h s as a function of the per- 
pendicular IC coupling Jq C (cf. Fig.jljD) and A = 1. (d) Phase 
diagram with critical IC coupling in 3D and 2D (thin line). 
Phase boundaries are extracted from the kinks in the h s as in 
(a-c). Red dashed lines: analytical HCB results [Eq. ( |S54[ )]. 
Symbols: the dependence of the critical Jq C on the uniaxial 
exchange anisotropy A — 1 in 3D for a = 0.5, where o/x cor- 
respond to the DMRG/analytical HCB results, respectively. 



Figure 4: Phase diagram for MPP swith skew (diagonal) IC 
coupling J\° (left) and j\° (right) in 3D. 



on A, a, and N\q. For the nematic phase we obtained 
exact values of h s using the HCB-approach|2"6"I |2"? . The 
HCB values are in full accord with the DMRG results. 
In the limit J IC « 1 we arrive at the analytical ex- 
pansion h S) 3 ~ hl D + ri 2 {J lc ) 2 + r7 4 (J IC ) 4 , which is ap- 
proximate but accurate enough for our present purposes 
and where h 1 ® = -A + 2a + A 2 /(2A + 2a) g6J, and 
2 m (a,A) = N IC (A+a)(3a 2 +3aA+A 2 )/[A{A+2a)} 2 w 
Nic (5/6 + 3a/4), when A ~ 1. The expression for 
the next, quartic term 774 is provided in Ref. 1271 Com- 
paring the expressions for h Sj i and h s ,2 one notices the 
presence of nonlinear IC terms and a two times smaller 
linear term in the nematic phase as compared to the 
usual 1-magnon phase. The solution of the equation 
h s ,2 = h St i gives analytical expressions for the critical 
IC interaction Jqcj- Keeping only the linear term in 
the expression for h s 2 , we find (cf. Eq. (51) in Ref. [TT)l 
|jic ri | = (4aA 2 -A-a)/[4a(A + a)iVic] and including 
the quadratic term \27\ , we obtain 
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A comparison of the numerical DMRG results in Fig. 
[3] (cf. Fig. 6 of Ref. |3S) shows that Eq. ( |S54"| is very 
accurate for 3D systems and works well for 2D ones, too. 

The phase diagram for the situation of the two other, 
skew types of IC interaction, j{ c and J^ c (see Fig. [it 
and d) are shown in Fig. |4j An inspection of the phase 
diagrams reveals that the maximal value for the critical 
J IC always occurs in the nematic phase at a slightly be- 
low 1, i.e. in the region of maximal in-chain frustration 
and quantum behavior |26l [2"5] . For the situation of per- 
pendicular coupling this can be understood already in 
linear approximation, where j cr ,i is proportional to the 
difference of 1- and 2-magnon critical fields of an isolated 
chain j CTi i = 2(/i^ D — /ig D 1 )/7Vic. Near the critical point 
(a > 1/4) and for almost decoupled Heisenberg chains 
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Figure 5: Phase diagram and pitch (contour lines) as a func- 
tion of the diagonal IC coupling Jl c (in units of |Ji|) and 
the uniaxial exchange anisotropy A — 1 for a — 0.36, as is 
relevant for linarite, PbCuS04(OH)2- The red contour line 
corresponds to the experimental value of the pitch, 34°. 

the saturation field tends to the simple 1-magnon value 
and additional quantum effects vanish. 

Having investigated theoretically in general how the 
competition between frustration, different types of IC 
coupling and exchange anisotropy plays out, we now 
apply these insights to identify candidate materials po- 
tentially displaying a quantum MPP. Li2Cu02 is near 
the critical point, having a ~ 0.33 and a rather small 
A-Ik 0.01 [24]. Its IC coupling j\ , however, is strong 
enough to even destabilize the spiral state and drives the 
chains FM. Also Li2ZrCu04 is close to the critical point 
(a ss 0.3 [17]) but in this case as well for any realistic IC 
interaction and reasonable value for A, all higher MPP 
are unstable. The compounds Li(Na)Cu202 are away 
from the detrimental critical point but their IC coupling 
is too large (J IC ~ 0.5 to 1 |4"0"H4"2"] ) to establish a ne- 
matic phase for the estimated, moderate, values of A |43| . 

Instead LiVCu04 is a good material for a nematic 
phase, having a coupling between the chains that is char- 
acterized by a very weak Jg C , which manifests itself in 
strong quantum fluctuations evidenced by a small or- 
dered magnetic moment (0.3/i B ) at low temperature and 
the observation of a 2-spinon continuum in inelastic neu- 
tron scattering [44] . The weak J IC is also in accord with 
the fact that its saturation field is close to the value of 
the uncoupled ID-chain given by /iJ D [45]. In addition, 
the estimated a « 0.75 [28, 45 , near the maximum of the 
critical Jo^ r (l/a)-curve is almost optimal for a nematic 
phase to survive (see Fig. [3| . 

A very interesting case is provided by the natural min- 
eral linarite, PbCuS04(OH)2, which consists of neutral 
edge-shared Cu(OH)2-chains surrounded by Pb 2+ and 
[S04]~ 2 ions and has a 0.36 [21] . Below 2.7 K a spiral 
state with a pitch of 34° sets in [221 [46] . A perpendic- 
ular Jq C barely affects the pitch of the spiral, in sharp 
contrast to skew IC coupling J\ G . We have considered 



this situation theoretically in more detail and calculated 
the phase diagram as a function of A — 1 and j{ c , see 
Fig. [H] For the given value of a a small j[ c and A — 1 
are enough to reduce the pitch from about 60° to the ex- 
perimental value of 34° . The experimental pitch strongly 
restricts the possible values for j\ c and A — 1 (see the 
red line in Fig. [5] ) . An additional piece of information 
is the experimental value of the saturation field of 1 1 T 
- the ID saturation field gives in this case about 5 T - 
which indicates a reduced value of J} c , renormalized by 
a sizable A — 1, placing the system close to the triatic, 
3-magnon region of the phase diagram in Fig. [5] 

In this context experimental studies under chemical 
or physical pressure are of great interest, since these can 
significantly change the IC coupling. When applying hy- 
drostatic pressure one expects an increase of the IC cou- 
pling and thereby a weakening and possibly disappear- 
ance of the MPPs in the mentioned two candidate mate- 
rials. Vice versa, growing isomorphic crystals with larger 
isovalent cations, i.e. substituting e.g. Li or Na by Na, 
Rb, or Cs, respectively, is expected to lead to candidate 
MPP materials due to a decrease of IC couplings. If pos- 
sible to synthesize one expects e.g. for Cs(Rb)Cu202 
and Na(Rb) 2 ZrCu04 an increased stability of the ne- 
matic and triatic phase, respectively. Preparing strained 
epitaxial thin films from candidate materials will cause 
similar effects, where a tuning of the strain can change 
the IC in different directions. 

We have, in summary, demonstrated the crucial role 
of different types of antiferromagnetic inter-chain interac- 
tions and the uniaxial exchange anisotropy in frustrated 
quasi-lD helimagnets. The rich and exotic physics of 
multipolar phases recently predicted for single chains is 
very sensitive to the strength and type and these ad- 
ditional and unavoidable interactions. Unfortunately, 
this prevents a realization of multipolar phases in most 
presently known spin-chain materials. But we find at 
least two notable exceptions: LiVCu04, where a nematic 
phase is expected, and linarite, PbCuS04(OH) 2 , which 
according to our present calculations is in the close vicin- 
ity of a triatic instability. In addition we proposed several 
new material systems as potential candidates with mag- 
netic multipolar ground states and point out the large 
experimental potential of tuning the interchain interac- 
tions by pressure and strain. 
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RI615/16-1 (JR)] for financial support and H. Rosner, 
A. Wolter, and M. Schaeper for discussions on linarite. 
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We provide details on the derivation of the equations in the main text, following the approach 
developed in Ref. 26. The calculations are tedious but straightforward. 



At high magnetic fields, the Hamiltonian of coupled frustrated spin-1/2 chains with the ferro- antiferromagnetic 
J\- J2 XXZ-Heisenberg model reads 



H — Hid + Hi, 



Hid — 



E 



2 

(MS, 



A QZ QZ 1 A+ &- 



Hie - ^ X/ ^ f ^m^m+f + ^tn^tn+i 



(SI) 
(S2) 

(S3) 
(S4) 



where m enumerates the lattice sites, Is r = ±na, n — 1, 2 determines the NN sites within the chain, and a is the 
lattice vector along the chain. The vector f connects sites at different chains. We restrict ourself to the case of uniaxial 
exchange anisotropy and the magnetic field directed along that axis, fi = g^ B . 
In terms of hard-core boson operators b, defined by 

S+ = b, §- = &t S z = 1 -h, 



n m = bl n b m = 0, 1, 



{^m,&L} = 1 , 



— 0, m ^ m', 

tmttt ' ' •) 



the Hamiltonian (SI) becomes 



bl\FM) = bl\---f 

= i---m m m 

(bl) 2 = (b m f=0, 



H — H + H int , 

m m,R 

Hint = ^ E ^R^R-^m^m+R: 
m,R 



(S5) 



(S6) 
(S7) 

(S8) 



where w = fM-~ J2n 4Ar, R = r, f 

The n-particle excitation spectra are given by the singularities of the corresponding retarded Green's functions 
(GF) 



((X\Y)) = -if dte 1 ^ 1 -^ 



u>((X\Y)) 



X,Y 



X,H 



\Y))- 



(S9) 
(S10) 



A negative value of the excitation energy signals an instability of the ground state, which is given by the fully polarized 
state achieved for a magnetic field above the saturation field T~L > H s - 
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Figure SI: Cartoon of the effective impurity problem given by the Hamiltonian ( S 13 1 , which describes the internal motion 



of a magnon pair with the total quasi-momentum k. The pink, open, shaded and cyan circles depict the impurities with 
£m = oo, J1A1, J2A2, J±A± respectively, • : the regular sites of the lattice, arcs: the k-dependent hoppings. 



The equation of motion for the two-magnon operator 
reads 



(2uj + J2r >/rAr(5i,r) A k ,i 

(1 - Sl,o) E R J R cos ^ ^k,l+R, 



(Sll) 



(S12) 



where k being the total quasi-momentum of the magnon pair, N = Nj_N x is the number of sites, N± is the number 
of chains, and N x denotes the number of sites in the chain. 

As usual, the exclusion of the center of mass motion reduces the problem of an interacting pair particles to a one- 
particle problem of motion in an effective potential well (EPW) . In our case it corresponds to an impurity problem 
in a tight-binding Hamiltonian [55] (see Fig. Sll 



H tb (k) = f(k) + V, 
f(k) = 2 Wo ^|m)(m| 



|m + R)< R (k) (R| 

m,R 



V 



E 



m e r 



(m 



(S13) 
(S14) 

(S15) 
(S16) 



where 



t R (k) = Jrcos 
m' = 0, r, f s 

The Hamiltonian depends on the total pair momentum. 
The two-magnon GF reads 



kR 

OO, £r 



JrA 



R^R- 



G 1)n (k,w) = ({A Kl \Al, 

= (<h\ (u) - H^y 1 \<t> n ) 



(S17) 
(S18) 



(S19) 
(S20) 



with \(f>i) = (|1) + |— 1)) /v2- The GF is analytic everywhere in the complex energy plane but may have singularities on 
the real axis: branch cuts and isolated poles. The branch cuts correspond to the continuum spectrum of unbounded 
motion of the effective particle, which in its turn correspond to the two-particle continuum in the pair motion. The 
poles correspond to the energies of localized impurity states, which are bound states for the pair when the energies 
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lie below the continuum or anti-bound states in the opposite case. It is clear from Eqs. ( S13 (-( S18 1 that bound 



states are possible only when some £r are negative, i.e. for FM Jr < 0. The bound state energy and the continuum 
boundaries depend on the total momentum of the pair k. If the bound state energy minimum lies below the lowest 
continuum energy (that may occur at different k- values), the bound pairs will condense at magnetic fields just below 
the saturation field, the gas of pairs being the nematic state of the magnetic svstem[51 \TT\ . 

When all Jr are positive, like in AFM-AFM J\-Ji model, only anti-bound states occur at energies higher the 
two-particle continuum. In this case only the one-magnon condensation occurs below the saturation field. 

We will use the identity 

G = g + gVG, (S21) 
for the solution in the real space of the impurity problem given by Eqs. ( S13 )-( S20 1 (see Fig. SI I. In Eq. ( S21| , 



g = (qj — Tj is the resolvent operator for the periodic part, and G = (uj — H t bJ is the resolvent for the impurity 
problem. According to Ref. 57] we may solve the problem step by step. Starting from the GF of a free particle, which 
in the matrix form reads 



9l,n = 51-n(k,Cj) 

1 >— ^ cosq(l — n) 



sir 



UJ 



i , ,sw i , ,sw 

- L " 1 W k/2+q + W k/2-q 



(S22) 
(S23) 

(S24) 



R 



we add the impurity at the origin. Its infinite potential reflects the impossibility to have two particles on the same 



site (S5) 



9\ 



(0) 



(°) 

.9i,n + .9i,oeoffo,ni 



(0) . Sl0 e 0fl0, n 9l,090,n 
.9l,„ = 9l,n + \ T^T > 9l,n 



1 — £o9o,o go,o 

Next, we add an impurity at the site i and express the GF via g^ 

(o) (o) 

(i) _ (0) 9\,j £ i9-un 
5l, n — 9 1,1! + " (0) ' 

1 - e i.9i,i 

and so on, the GF of the system with r impurities is expressed via the GF of the system with j — 1 impurities 



(S25) 



0) 

9l,n 



(r-1) 



(r-1) 
.9l, r £r.9r,n 

(r-1) 
1 - £r.9r,r 



(S26) 



Thus, in principle, we may take into account any number of in-chain and inter-chain exchange couplings (IC) and 



obtain Gi jn (k,w) ( S19 1 . The explicit expression for the GF G± t i(k,u>) for the ID J\-Ji model (S3) has been given in 



Ref. 1261 It's spectral density is plotted in Fig. S2 The sharp k-dependent peaks below the two-particle continuum 
corresponds to bound pairs of magnons. 



At higher dime nsions, the role of the inter-chain interaction (|S4|) is t wofold . First, the periodic part of th e effective 

Second, 



Hamiltonian (|S14|) becomes D-dimensional. This changes g from Eq. (|S22| via the change of oj^. w (IS24 1 



new impurities with the strength e r = Jj_Aj_ are added at points r. The simplest geometry for the IC corresponds to 
f-vectors perpendicular to the chains, which connect NN sites, only. The spectral density for GF G a , a (k, ui) for k || a 
for the 2D case is depicted in Fig. [S3] We see that for small IC couplings the spectral density behaves qualitatively 



similar to the ID case, i.e. the peak corresponding to the bound pair lies below the continuum (left panel of Fig. S3 1, 
and its dispersion exhibits a minimum at the total momentum ka = tt of a pair. We have checked numerically that 
the minimum position remains at the point k,,- = (vr/o, 0, 0) for all values of IC satisfying the condition J± < J cr . On 
the right panel of Fig. [S3] we see that the behavior of the spectral density changes for large enough IC. The bound 
state is still present near the edge of the Brillouin zone, but its energy is higher than the minimum of the two-particle 
continuum. It is clear that the critical IC value J cr is defined by the condition 



LO h = ^(kjr) = W v 



(S27) 
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G>2uH 

Figure S2: The spectral density of the two-particle Green's function for an isolated chain. ID case, i.e. Ji=-1, J2=l, Jj_=0. 
Cyan and magenta thin lines shows the lower boundary of the 2-magnon continuum. 



J 1= -1, J 2 =1, J y =0.1 




-2.8 -2.6 -2.4 -2.2 -2 -1.8 -1.6 -1.4 -1.2 -1 
co-2nH 



25 
20 
15 
10 
5 




J 1= -1, J 2 =1, J y =0.25 



IrnfO, j(k,co-2uH)]/7i+a*(7i-k) 




Figure S3: The spectral density of the two-particle Green's function for a 2D arrangement of unshifted J\-Ji chains and 
perpendicular IC interaction, left: J y = J± < J cr , right: J y > J cr . 



where oj r . 



2 (fiH — \ Ji\h Si i) is the minimum of the energy of the two-particle continuum, and 



h Sl i = -Ai + a(A 2 + l) + 

+ 0.125/a + 0.5AT io (j ic A ic + |j lo |) , 



(S28) 



is the critical field of the 1-magnon instability (Eq. (1) of the main text). In order to find the expression for the 
saturation field T-L s as a function of IC |Jj_| < J cr , we need the expression for Wf,, which is the position of an isolated 
pole of the GF 



[Ga,a(k*,W6)] '=0. 



(S29) 



In terms of the effective model //^(k^) (S13), u>b is the energy of the localized impurity level. From Eq. (S17) we see 
that the nearest-neighbor hopping along the chain vanishes t a = J% cos J = 0, and the sites with r = na + mb + Ic 
having odd and even ra's are decoupled. In the subsystem with odd n's, only two impurities of the same strength 
e a = Ji are present at the sites ±a = (±a, 0, 0). The effective particle motion is not affected neither by the impurity at 
the origin (of infinite strength) nor by the impurities at the sites f = (0, ±6, 0), (0, 0, ±c) with the energies J2A2, and 
Jj_A_i_, respectively. Note that this peculiarity has an important consequence: the critical value of the IC given below 
by Eqs.(S54)-(S57l depends only on the nearest-neighbor exchange anisotropy value Ai. So, we may immediately 



write down the expression for the GF (cf. Eq. (49) of Ref. 



G a , a (k w , W ) = \g£1(K,u) 



- JiAj 



(S30) 
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where 



G^l(k n ,u,) = (0 a |( W -T(kJ-|O) £o (O|) Va) 



= goikn) + g 2a (k w ) 

= 5o(kir) + 32a (K)- 



9o 



(S31) 
(S32) 



In Eq. (S31) we have used the relation ( S25 1 and Eq. ( S32 1 follows from g a (k T ) = 0, since the vector a joins two 
decoupled subsystems. Then Eq. ( S29 1 may be rewritten as 



Now, using the definition ( S22 1 , we may write 



N 



1 >— ^ 1 + cos 2q x a 

uj-E 1D (ir,q x ) [ 



Ei D (ir,q x ) = 2 [yM + Ji (cos q x a - Ai) 
+ J 2 (cos2q x a - A 2 )] , 
E ± (n,q) = N tc J ± ( 7q - A ± ) , 



(S33) 

(S34) 
(S35) 



(S36) 
(S37) 



where 7 g = cosq y b ((cos q y b + cos q z c) /2), N[ c = 2(4) for a 2D (3D) geometry, respectively. In the 2D case the 
summation over q z should be dropped. The ID GF as given by Eq. ( |S35[ ) is easily calculated 



G${n,u) = G(z)/J 2 , 

G(z) = [z + l-r^f 1 , 

where we have introduced the dimensionless variable 

z(u) = [w-2{iM- JiAi - J 2 A 2 )] /J 2) 



(S38) 
(S39) 

(S40) 



and the dimensionless Green's function of a semi-infinite tight-binding chain t(z) = [z — t(z)] 1 . Now, we search for 
the solution of Eq. ( S30 1 in the form 



= J-2 1 -/,! + O + 2 ( iM - JiAi - J 2 A 2 - iiV ic JxAx ) 



where £ is unknown, and 



z bl 



Ai + a a 

— + x 

a Ai + a 



(S41) 



(S42) 



is the solution for the ID-problem 
a = J 2 /| Ji| as compared to Ref. ,261 



Note that here we use another definition for the frustration parameter 



Assuming ( <C 1, we rewrite Eq. ( S33 1 in the form 

1 



E\ -> lj' (/- — \ m — 

q y ,q z m=0 



where e q = N ic Jj_jq/J2, 
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Figure S4: (Color online) The saturation field h Bl 2{o-,j y ) = fM s /\Ji\ for a 2D array of chains (J y = J±) for a = l,Ai = 1. 
Black solid line: the result of analytic Eq.(S52l, green short-dashed line: the result of the expansion (S52| up to second order 



(i.e. £4 is neglected), black da shed line: the field of the 1-magnon instability h s ,i ( S28 1 . Points: DMRG-data and data from 
the numerical solution of Eq. (S29l. 



Note that G(zn) = — a/Af , and keeping only terms with m < 4, we obtain the equation 



2 1 el ) G" 



^(C 3 + 3Ce|J G'" 



24 



(S43) 



where 



— V e m 



N ± 



and we have taken into account that e q = e q = 0. The direct calculation yields e q = ( J±/ J2) 2 , and e q = 6 (J1/J2) 4 
(36(Jj_/J 2 ) 4 ) for 2D(3D) respectively; 

G' = G 2 [t'-1] 

G" 

G'" 



G IV = 



2G 3 y - if + G 2 t", 
6G 4 [t' - l] 3 
6G 3 [r' - 1] t" + G 2 r'", 
24G 5 [r'-l] 4 + 36G 4 [r'-l] 2 r" 



6G 3 (r") 2 + 8G 3 [t 1 - 1] r'" + G 2 t iv , 
a 2 



(S44) 
(S45) 

(S46) 
(S47) 



Ai (Ai + 2a)' 



-IV 



= -2 
= -6 
= -24 



a (Ai + a) 
Ai (Ai +2a)_ 

a(Ai + a) I 4 A? + 2Aia + 2a 2 



Ai (Ai + 2a) 
a 5 (Ax + a) 5 F 



Ai (Ai +2a) 



Substituting the expansion 



[Ax (A 1+ 2a)] 7 ' 
F = A 4 +4A 3 a + 9A 2 a 2 + lOAiCv 3 + 5a 4 . 

C = (die + (.ijic + (ait + (dL 



(S48) 
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(Jic = J±/\Jt\) into ( S43 1 , we obtain (i = £3 = 0, and 

N ir G" 



<2 



C4 



2a 2 G' ' 
N ic (Ai + a) 

a[Ai (Ai + 2a)f 
_JL_ l(T 
~~G' ~2~ 



[A 2 + 3Aia + 3a 2 ] , 



ci 



N ic G" 
2a 2 



<2 + 



G IV 

24/3 4 



At the saturation field, the Wb in the right-hand side of Eq. (S41 ) vanishes, and we obtain 



h s ,2 = h 



ID 

s,2 



: ii C A^ 



(CaiS + C4Ji 4 c ) 



(S49) 
(S50) 

(S51) 
(S52) 



.id 



-Ai + oA 2 - -z b i 



where /i s = fjfH s /\Ji\. Eq. (S52l coincides with Eq. (3) of the main text with rji = —a(i/2. Its validity is demonstrated 
in Fig. |S4| As an example, we have chosen the 2D case and a = 1, i.e. the optimal region for the existence of the 
nematic phase, where j cr w 0.167. For small j± = J±/\Ji\ the second order expansion reproduces well the DMRG 




Figure S6: The saturation field h S: 2 for the 3D case for a = 0.5 (J y 
taken into account. The mean ing of the lines is the same as in Fig. 
numerical solution of Eq. ( S29 1 (Ai = 1) 



= J_l). The easy-axis anisotropy of the NN coupling is 
Si] Points: DMRG-data (Ai / 1), and the data from a 
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data which coincide with the results from a numerical solution of Eq. ( S29 ) . Naturally, for larger interchain coupling 
j± > 0.15 the fourth order expansion is needed. 

The boundary between the 1-magnon and the 2-magnon phases is obtained by solving the equation h s ^(jci) = 



h s i(j CI ) for the critical IC j cr . If one retains only the linear term in the expansion in powers of the IC given by (S52 1 
we obtain (cf. Eq. (51) in Ref . [TT)| 



licr,l| 



AaAj - Ai - a 
W ic a (Ai + a) ' 



(S53) 



This approximation demonstrates the qualitative behaviour of j CI as a function of the anisotropy and the frustration 
parameters Ai and a, respectively. Practically, a fully quantitative agreement with our numerical data is achieved, if 



we account also for the quadratic term in Eq. (|S52|) 

\jcr.: ' 



2a( 2 



N lc + ^Nl + 4N lc a( 2 \jcr,i\ 



(S54) 



It is convenient to normalize the couplings on J 2 > 0, and introduce n = 1/a, which measures the attraction provided 



by the FM J\. Using the same normalization for the IC, too, we write y = J±/J 2 = jic/oi. Then the Eqs. ( S50 1 

N lc (/tAi + 1) 



( S53 1 , and ( S54 1 may be rewritten as 

C2 



[Ai (kAx + 2)]' 



[k2AI + 3Aik + 3] , 



\Vcr,l\ = 



\Vcr.2\ 



2 (4A? - kAi - 1) 
AN ic («Ai + 1) 



-Ni, 



A :2 



4iV lc C 2 |y, 



cr.l | 



(S55) 
(S56) 
(S57) 



A comparison of the results of the approximate analytic Eqs. ( S56 1 and ( S57 1 with the numerical data is shown in 



Fig. S5 Note the high accuracy achieved already in the second order of the IC in Eq. (S52). Finally, an example of 



the saturation field dependence on the anisotropy parameter is shown in Fig. |S6| 



